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We study the coarsening dynamics of two and three dimensional biaxial nematic liquid crystals, 
using Langevin dynamics. Unlike previous work, we use a model with no a priori relationship 
among the three elastic constants associated with director deformations. We find a rich variety 
of coarsening behavior, including the simulataneous decay of nearly equal populations of the three 
classes of half-integer disclination lines. The behavior we observed can be understood on the basis 
of the relative values of the elastic constants and the resulting decay channels of the defects. 

PACS numbers: 64.70Md, 61.30Jf 

Topological defects play an important role in the equilibration process following a quench from a disordered to an 
ordered phase ("coarsening dynamics"). Coarsening dynamics in nematic liquid crystals, particularly with uniaxial 
ordering, has been the subject of much active investigation in recent years in theory, experiments and simulations 
[1, 2, 3, 4], in part because of the rich defect structure of liquid crystals. On the other hand, relatively little attention 
(with the exception of the two-dimensional work of Ref. [2]) has been paid to coarsening dynamics in biaxial liquid 
crystals, in part because of the dearth of experimental realizations of biaxial liquid crystalline phases. However, biaxial 
nematics have many unusual topological features, which might be expected to influence their coarsening dynamics and 
thus warrant study. Biaxial nematics differ from their uniaxial counterparts in that they possess four topologically 
distinct classes of line defects (disclinations) , while possessing no stable point defects (except in two dimensions where 
the line defects reduce to points) [5, 6]. The classes of disclination lines are distinguished by the rotation of the long 
and short axes of the rectangular building blocks of the system. In the first three classes one of the three axes is 
uniformly ordered, while the remaining two axes rotate by 180° about the core of the defect. The fourth class consists 
of 360° rotations of two of the three axes. The disclination lines form closed loops in three dimensions (with a single 
defect class per loop) or form a network where three lines, each from a different class meet at junction points [7]. The 
fundamental homotopy group of biaxial nematics is non-Abelian leading to a number of interesting consequences. 
E.g., the merging of two defects will depend on the path they follow, and two 180° disclinations of different types will 
be connected by a 360° "umbilical" cord after crossing each other [8]. 

Zapotocky et al. [2] studied coarsening dynamics in a two-dimensional model of biaxial nematics, utilizing a cell 
dynamical scheme applied to a Landau-Ginzburg model, where the gradient portion of the energy was given by: 
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Here M is a coupling constant, and Q a p is the symmetric-traceless nematic order parameter tensor. Repeated 
indices are summed over; in the case of a two-dimensional nematic, a is summed over x and y, while [3 and 7 are 
summed over x, y and z. Zapotocky et al. found that of the four topologically distinct classes of disclinations, only 
two classes (both corresponding to "half- integer" lines, i.e., 180° rotations) were present in large numbers at late 
times. Subsequently, Kobdaj and Thomas [9] showed within this one-elastic constant approximation that one class 
of half-integer disclinations is always energetically unstable towards dissociation into disclinations of the other two 
half-integer classes. 

In this Letter we show that if one considers a more general gradient energy the coarsening dynamics of biaxial 
nematics is much richer than what occurs with the above simple model. In particular, with appropriate sets of 
parameters one can obtain a coarsening sequence with all three classes of half-integer disclinations present in nearly 
equal numbers even at late times, or a sequence with only one class of half-integer disclinations surviving until late 
times. When all three classes are present, the topology of the coarsening sequence in three dimensions is markedly 
different from the uniaxial case. 

To understand why the model free energy of Eq. (1) is not general enough for coarsening studies, it is helpful to 
see what it yields for the director elastic constants. For biaxial nematics there are three directors which form an 
orthonormal triad of vectors u, v, w describing the alignment of the constituent "brick-like" molecules. The tensor 
Qap can be written in terms of the orthonormal triad as: 

Qaf3 = S(w a Wj3 - -5 af 3) + T(u a Up - V a Vp), (2) 

where S and T are respectively the uniaxial and biaxial order parameters. If we insert Eq. (2) into Eq. (1), we find 
[10]: 



Fgrad = -M(d a Q^)(d a Qp 7 ). (1) 



2 



F glad = K u [(u ■ Vv • w) 2 + (v Vv • w) 2 + (w Vw • v) 2 ] + K v [(v • Vw • u) 2 + (u • Vu • w) 2 + (w Vw • u) 2 ] 

+ K w [(yr- Vu-v) 2 + (u- Vu • v) 2 + (v Vv u) 2 ], (3) 

where 

K u = 2M(S + T) 2 
K v = 2M{S-Tf 

K w = 8MT 2 (4) 

Thus, there are three elastic constants K u , K v , K w , each corresponding to one of the three classes of line defects. Each 
class corresponds to rotations of two of the three vector fields u, v, w, about the defect core, with the third vector of 
the triad undistorted. We denote the classes as follows [6, 7]: C u (u undistorted), C w (w undistorted), and C v (v 
undistortcd). The energy of a defect in class C'i, i ~ u, v, w is proportional to the elastic constant K{. Note, however, 
that in the model specified by Eq. (1) only two of the three elastic constants are independent, as they are related as 
indicated in Eq. (4). In fact, the specific relationship among the three elastic constants gives rise, as shown in Ref. 
[9], to the presence of only two defect classes at late times. Irrespective of the values of S and T (M simply sets the 
overall scale of all three elastic constants), one of the three elastic constants is always greater than the other two, 
yielding a decay channel for the defect in the class with the largest elastic constant. 

There is no symmetry reason to restrict our attention to the model free energy, Eq. (1). Even if we neglect the 
clastic anisotropy associated with bend, splay and twist distortions, a biaxial nematic should be described in general 
by three independent elastic constants, K u , K v and K w [10]. In terms of the order parameter tensor Q a p, this requires 
a term of third order in i^rad : 

-Fgrad = ^M{d a Q M ){d a Qf }l ) + LQ a p(d p Q ai )(d p Q 01 ), (5) 

which upon substituting Eq. (2) yields the elastic constants: 

K u = 2(M(S + T) + LST)(S + T) 
K v = 2(M(S-T)-LST)(S-T) 

K w = 4{2M-LT)T 2 (6) 

Because of the extra coupling constant L, there is no predetermined hierarchy among these elastic constants. Unlike 
the model of Eq. (1) it is now possible with appropriate choices of S, T and the ratio L/M to have all three elastic 
constants equal (which will give rise as we demonstrate below to a coarsening sequence with nearly equal populations 
of the three classes of half-integer defects), or have one constant smaller than the remaining two, yielding a coarsening 
sequence with only one class of defects at late times, or recover the behavior seen in Ref. [2]. 

To simulate the coarsening dynamics of this general model of biaxial nematics, we consider its lattice analog 
introduced by Straley [11]. In this model the interaction between two biaxial objects located at sites i and j of a 
cubic lattice with orientations specified by the orthonormal triad u, v, w is given by: 

Ua = ~0{wi ■ w,) 2 - 2 7 [(u 4 • u,) 2 - ( Vi • v,) 2 ] - 

6 -[(u t ■ u,) 2 + (v 4 • v,) 2 - (u, ■ v,) 2 - ( Vi • u,) 2 ]. (7) 

This model has a phase diagram with two uniaxial phases, one with rodlike order (alignment of the w vector field), 
one with discotic order (alignment of the u vector field) and a biaxial phase with alignment of all three vector fields 
[11, 12, 13]. 

The elastic constants emerging from Eq. (7) can be determined by considering the interaction between two objects 
which arc aligned in turn along each of the three directions u, v, w, with the results: 

3/3 5 

K w = 25 (8) 
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As in the continuum model Eq. (3), the parameters 0, 7 and 8 give rise to three independent clastic constants. Stability 
requires that 3/3 + S > 4\-f\, and 0,8 > 0. 

We have simulated the coarsening dynamics associated with Eq. (7) using Langevin dynamics, expressing the three 
unit vectors w, u and v in terms of Euler angles cj>, 9 and ip. The equations of motion are given by : 

where U is the sum of Uij over the nearest-neighbors of site i, £ is a damping coefficient and Rj>(t), Rg(t) and R#(t) 
are uncorrelated random thermal noise sources. Each of the random variables R has a Gaussian distribution of 
variance 2ksT(/St where ks, T and St are Boltzmann's constant, the temperature and the time step respectively. We 
measured time in units of (choosing ( = 1), and used a timestep St — 0.0005. A dynamical equation for cos 9 rather 
than 9 must be used in order to reach the correct equilibrium states [14]. We verified that our dynamical equations 
led to the same phase diagram produced by Monte Carlo simulations [12, 13]. 

As in previous numerical studies of defect behavior [15, 16], we introduce a disclination line segment counting 
operator, 

DTjkl = \\ l - sgn{(w 2 • Wj)(wj • w fc )(w fc • Wi)(wj • wj)}] (10) 

which is unity if a disclination line segment pierces the lattice square defined by the four vectors Wj , Wj , wj, and w; . 
We define analogous operators D^- kl and D^- kl for the u and v vectors on this lattice square. In principle either two 
or none of the three operators should be unity for a given square, and thus we can assign the line segment to one 
of the classes C u , C v or C w . In practice, we found a small number of squares where only one or all three operators 
were unity, an artifact of the discreteness of the underlying lattice. We obtained physically reasonable results by 
classifying the defects on the basis of the operators D^- kl and D^ kl , assuming that DJ- kl was unity only if one of the 
two former operators was unity. This procedure always yielded closed defect loops in three dimensions, a reasonable 
test of our algorithm. We determined the location of the integer-valued defects where either w or u (or both) rotate 
by ±360° degrees using the algorithm of Ref. [4] which provides an upper bound on these defects. We found very few 
integer-valued defects using this method, so a more accurate algorithm is not needed. 

We quenched the system instantaneously from an initial configuration of random orientations of the three vectors 
u,v,w (i.e., a high temperature state) to zero temperature (a biaxially ordered state for all S ^ 0), and then let the 
system evolve in time according to the dynamical equations (9) monitoring the defect populations both visually and 
statistically. We studied the model in both two and three dimensions. 

In agreement with our arguments above regarding the elastic constants, we found three qualitatively distinct types 
of coarsening behavior depending on the values of 7 and 8 (we set = 1). For S <C 1, and 7 satisfying the stability 
requirement given after Eq. 8 defects of class C w are energetically favorable compared to those of the other classes, 
and thus only these former defects survive until late times. For larger values of 5, but still less than 27 (with 7 > 0) 
we found defects belonging only to classes C u and C w , consistent with Eqs. (8), namely, K v > K u + K w , and the 
results of Ref. [2]. This regime also includes the parameterization of Eq. (2) used in Rcfs. [12, 13], where 5 = 47 2 /3 
was chosen. In Fig. 1 we show the total line length of defects in each of three classes as a function of time after the 
quench. In Fig. 2 we show a snapshot of the simulation cell late in the coarsening sequence. The C v class defects 
have disappeared and the C u and C w class defects remain as nonintersccting closed loops. When 7 < and S < 2\"f\, 
we have K u > K v + K w , and we found, consistent with this inequality, only defects of classes C v and C w at late 
times. We found no apparent line crossings, entanglements or junction points (the latter obviously would require 
three classes of defects). The two classes of loops appear to coarsen independently, as can be seen in the animations 
on our web site [17]. 

The most interesting and novel coarsening sequence occurs when 7 = and 8 = 1. For these parameters Eq. (8) 
implies that all three elastic constants are equal, and we found a coarsening sequence with nearly equal populations 
of the three defects. This behavior persists over a range of parameters S ~ 1 ± 0.5 and 7 ~ ± 0.2. The coarsening 
sequence for these parameter values is particularly interesting in three dimensions. Soon after the quench a uniform 
network of junction points where three disclination lines (one from each of the three half-integer classes) meet is 
formed. These junction points are illustrated in Fig. 3 at late times after the quench. The junction points are 
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distributed in a nearly uniform fashion throughout the simulation cell, and the distance between neighboring points 
grows on average with time [17]. The total line length of each of the three classes are nearly equal throughout the 
coarsening process (see Fig. 4). When the distance between neighboring junction points becomes comparable with the 
size of the simulation cell (see Fig. 3), the coarsening process is impeded. The final annihilation of the disclination lines 
can only occur via the shrinkage of individual loops. The formation of loops requires that some pairs of neighboring 
junction points approach each other, shrinking the line joining them while possibly increasing the length of the other 
two disclination lines attached to the pair of junction points. Ultimately, the pair of junction points meet at a "pinch 
point", where four disclination line segments corresponding to two defect classes meet. Subsequently the four line 
segments dissociate into two nonintersecting single class line segments as can be seen in our animations [17]. When 
this process has occured a sufficient number of times, individual disclination loops are formed which then shrink 
independently. 

In conclusion, we have shown that the coarsening dynamics of biaxial nematics is very rich, with late time behavior 
governed by either one, two or three classes of half-integer line defects, depending on the parameters of the system. 
When all three classes are present, a network of junction points is formed and a novel coarsening sequence occurs. 

We thank S. C. Ying, J. M. Kosterlitz and M. Zapotocky for helpful discussions, and G. B. Loriot for computational 
assistance. This work was supported by the National Science Foundation under grant DMR-9873849. Computational 
work in support of this research was performed at Brown University's Theoretical Physics Computing Facility. 
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FIG. 1: Log- log plot of the total line length of C v (o), C w (box), and C u (A) class defects for system size 40 3 , and parameters 
7 = 1/2, 8 = 47 2 /3 = 1/3, as a function of time after the quench. The data has been averaged over 60 initial random 
configurations. 




FIG. 2: The configuration of C„ (thin line) and C w (bold line) class defects for system size 40 3 , and parameters 7 = 1/2, 
5 — 4j 2 /3 — 1/3, at dimensionless time 80 after the quench. Note that all C v class defects have vanished at this late time, 
therefore C m and C u form nonintersecting loops which coarsen independently. 
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FIG. 3: The configuration of C u (thin line), C v (dotted line), and C w (bold line), class defects for system size 40 3 , and 
parameters 7 = 0, 8 = 1 , at dimensionless time 60 after the quench. 
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FIG. 4: Log- log plot of the total line length of C v (o), C w (box), and C u (A) class defects for system size 60 3 , and parameters 
7 = 0, 8 = 1, as a function of dimensionless time after the quench. The data has been averaged over 40 initial random 
configurations. While in principle the population of the three classes of defects should be equal immediately after the quench, 
there is a small difference in these populations due to the numerical issues discussed in the text following Eq. (10). 



